Numerical method for non-linear steady-state transport in one-dimensional correlated 

conductors 
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We present a method for investigating the steady-state transport properties of one-dimensional 
correlated quantum systems. Using a procedure based on our analysis of finite-size effects in a 
related classical model (LC line) we show that stationary currents can be obtained from transient 
currents in finite systems driven out of equilibrium. The non-equilibrium dynamics of correlated 
quantum systems is calculated using the time-evolving block decimation method. To demonstrate 
our method we determine the full I-V characteristic of the spinless fermion model with nearest- 
neighbour hopping tu and interaction Vh using two different setups to generate currents (turning 
on/off a potential bias). Our numerical results agree with exact results for non-interacting fermions 
(Vh = 0). For interacting fermions we find that in the linear regime eV <C 4£h the current I is 
independent from the setup and our numerical data agree with the predictions of the Luttinger 
liquid theory combined with the Bethe Ansatz solution. For larger potentials V the steady-state 
current depends on the current-generating setup and as V increases we find a negative differential 
conductance with one setup while the currents saturate at finite values in the other one. Both effects 
are due to finite renormalized bandwidths. 

PACS numbers: 71.10.Fd, 71.10.Pm, 71.27.+a 



I. INTRODUCTION 

Much of our understanding of electronic transport 
in solids is based on a picture of weakly-interacting 
charge carriers such as Fermi liquid quasi-particles. One- 
dimensional electron systems are well-known examples 
where this approximation fails. The low-energy physics 
of these systems is described by the Tomonaga-Luttinger 
liquid (TLL) theory and is realized in quantum wires 
such as carbon nanotubes. nanowires in semiconduc- 
tor hetcrostructures, or metal atomic chains on surface 
substrates*^— 

The transport properties of one-dimensional correlated 
systems have been extensively studied during the last two 
decades^— A major goal is to determine and understand 
the current-voltage characteristic / — V of various sys- 
tems made of quantum dots and wires. Most of these 
studies have been restricted to a regime where the en- 
ergy scale of the current-generating electromagnetic field 
or potential bias is small in comparison to the energy 
scale of the unperturbated systems (i.e., the band width 
in lattice models and the Fermi velocity in continuum 
models). The non-linear regime has mostly been inves- 
tigated in quantum contact problems, where the inter- 
action is confined to a small region of the system. For 
instance, non-linear current-voltage characteristics have 
been calculated within the TLL theory such as the power- 
law / ~ V a behaviour for the transport through a weak 
link££ or the exact I(V) curve for the current through 
a point contact in a fractional quantum Hall edge state 
device^^ However, the TLL theory is limited to low- 
energy excitations with linear dispersion and thus to po- 
tential biases V which are weak compared to the band 
width. Only recently the implications of a non-linear 
dispersion have started to be considered^ Thus current- 



voltage characteristics from the TLL theory are actually 
limited to a weak-bias regime. 

There are few works presenting full current-voltage 
characteristics with a voltage up to the largest energy 
scale of the system (for instance, see Refs.0, Till . and[l4) 
and none is concerned with one-dimensional correlated 
conductors. Thus transport properties beyond linear re- 
sponse are poorly understood. We believe that it is im- 
portant to attain a better knowledge of the non-linear 
transport properties in one-dimensional correlated con- 
ductors. First, the validity of weak-bias approaches can 
be confirmed only if one obtains some quantitative esti- 
mates of non-linear effects. Moreover, non-linear devices 
play a significant role in electronics and studies of non- 
linear dynamics are required to reveal the full potential 
functionality of quantum wires as electronic circuit com- 
ponents. 

In this work we develop and apply a method for in- 
vestigating the zero-temperature DC transport proper- 
ties of one-dimensional correlated conductors for poten- 
tial biases up to the order of the band width. For this 
purpose we study a well-known one-dimensional lattice 
model described by the half-filled spinless fermion Hamil- 
tonian with nearest-neighbour repulsion. Even though 
this model is exactly solvable by the Bethe Ansatz and 
the low-energy physics are determined by the generic 
TLL phenomenology^, its transport properties in the 
non-linear regime cannot be obtained analytically. 

To determine the transport properties we simulate the 
quantum dynamics of single chains which are driven out 
of equilibrium by a potential bias between the left- and 
right-hand halves of the chain and calculate the result- 
ing currents through the middle bond. The simulation 
of out-of-equilibrium quantum many-particle systems is 
one of the major challenges in computational physics. 
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Recently, a family of numerical methods has been de- 
veloped to simulate the real-time evolution of quantum 
lattice systems such as one-dimensional Hamiltonians 
with short-range interactions*^— The most prominent 
ones are the time-dependent Density-Matrix Renormal- 
ization Group (td-DMRG) and the time-evolving block 
decimation (TEBD) mcthodi 18 ' 19 Various flavours of td- 
DMRG have been successfully applied in studies of quan- 
tum many-body dynamics. In particular, they have 
proven to be promising tools for investigating electronic 
transport in strongly correlated nanostructures and one- 
dimensional conductors! 12 ' 13 ' 20 " — 

Surprisingly, the original TEBD metho d 18 ' 19 has not 
been applied to electronic transport problems yet. Both 
td-DMRG and TEBD methods can be described within 
a common mathematical framework, the matrix prod- 
uct quantum states . 17 ' 31 ' 32 Their accuracy and efficiency 
depend essentially on the amount of entanglement in 
the quantum system and should be similar. Admit- 
tedly, the TEBD algorithm is restricted to a small fam- 
ily of systems (one-dimensional Hamiltonians or lad- 
der systems with nearest-neighbour interactions only) 
while td-DMRG techniques are more versatile. How- 
ever, the TEBD algorithm is naturally parallelizable 
and thus fully scalable, which is a required feature in 
high-performance computing, while the efficient paral- 
lelization of the DMRG algorithms remains an open 
challenge! 33 ' 34 In this work we employ the TEBD method 
for computing the non-equilibrium dynamics of the spin- 
less fermion model. 

With the TEBD method we can compute the non equi- 
librium quantum dynamics of lattice models with a finite 
number of sites N over a finite period of time t. As we 
are primarily interested in determining the DC transport 
properties, TEBD results must be extrapolated to the 
thermodynamic limit N — > oo and to the steady-state 
limit t — > oo. However, finite-size and finite-time effects 
are very complex in these out-of-equilibrium quantum 
systems (for instance, see Ref.l35l) and extrapolations are 
difficult because we know very little about the scaling of 
currents with system size N and time t. For this rea- 
son we have investigated this scaling in an exactly solv- 
able classical model, the so-called LC line. Using this 
information we have developed a method for quantum 
systems which allows us to obtain reliable quantitative 
results for stationary currents I from numerical data for 
rather small system sizes and short simulation times. 

In this paper we show that our extrapolation approach 
allows one to determine the full I—V curves of interacting 
one-dimensional conductors using the spinless fermion 
model for illustration. In the non-interacting case this 
model can be solved exactly using the equation of mo- 
tion method. The outcomes of this special case confirm 
our extrapolation method and reveal a negligible numeri- 
cal error for the TEBD simulation results. For interacting 
fcrmions comparisons with predictions from the TLL the- 
ory and the Bethe-Ansatz solution confirm the validity 
of our method in the linear response regime. While for 



the linear regime the specific setup does not matter, it is 
highly decisive for the non-linear current-voltage charac- 
teristics. We basically distinguish between two different 
ways of creating a current flowiS In the first setup (I) we 
apply an initial voltage and calculate the ground state 
which has different particle numbers in its two halves 
and then let the system evolve with an overall equal on- 
site potential. In the second setup (II) we calculate the 
ground state without a potential difference but turn it on 
for the real time evolution. We have found that while the 
general shape of the current-curve as a function of time is 
dominated by classical effects (i.e., which are also found 
in the LC line model), the current- voltage characteristics 
are primarily determined by the chosen setup. For setup 

(I) the system shows a positive differential conductance 
for the full voltage-range and saturates at a finite value 
for very large potential differences. For the second setup 

(II) the linear response coincides with the one for setup 
(I) but for higher voltages we observe a negative differ- 
ential conductance. Both effects are also present in the 
non-interacting case and come from the finite bandwidth 
and the non- linear dispersion (i.e., the energy-dependent 
density of states) of the excitations ! 35 ' 36 

Our paper is organized as follows: In the next section 
we introduce the spinless fermion model, and the meth- 
ods used in this work. In the third section we investigate 
the classical LC line to explain the basic behaviour of 
the current in finite and infinite one-dimensional chains 
and derive our method to extrapolate the stationary cur- 
rent from finite system results. In the fourth section we 
describe finite-size effects and the convergence to a sta- 
tionary current for infinite system sizes, while the I — V 
characteristic of the spinless fermion model and compar- 
isons with exact results are shown in the fifth section. 
Finally, we summarize our findings in the last section. 
Some calculations are detailed in appendices. 

II. MODEL AND METHODS 
A. Spinless fermion model 

We consider a one-dimensional lattice model represent- 
ing correlated conductors driven out of equilibrium by a 
potential bias. For spinless fermions the Hamiltonian 
without a potential bias is 

N-l 

Ho =-t H Yl + 

+Vkg(^-|) (1) 

where in denotes the hopping amplitude between 
nearest-neighbour sites, Vh is the Coulomb repulsion be- 
tween spinless fermions on nearest-neighbour sites, and 
n j = CjCj. At half filling (N/2 fermions in the A-site lat- 
tice) this Hamiltonian describes an ideal conductor for 
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FIG. 1. One-dimensional conductor consisting of two coupled 
leads. A bias between its right- and left-hand halves is applied 
and the current is measured in the middle of the system at 
the junction between both sides (dashed connection). 



— 2tn < Vh < 2tn- It can be interpreted as a system of 
spin-polarized electrons. 

The low-energy behaviour of this lattice model is de- 
scribed by the TLL theoryj^^ The generic properties of 
a TLL are determined by two parameters: The velocity 
of elementary excitations (rcnormalizcd Fermi velocity) 
v and a dimensionless parameter K. From the Bethe 
Ansatz solution we know the relation between TLL pa- 
rameters and the parameters in and Vh of the spinless 
fermion model at half filling 
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where ol is the lattice constant. To drive the system out 
of equilibrium we a use step-like potential bias between 
the left- and right-hand halves of the chain 



(4) 



The potential energy step is set by Ae = |eV| where V 
is the voltage bias and e is the elementary charge. It 
is possible to use a smoother potential profile but the 
results for the stationary current using our extrapolation 
method are only slightly affected by the specific shape 
as long as the non-constant part in the middle is rather 
smooth and locally confined. 

The above system can be seen as two coupled interact- 
ing leads made of the sites {1 . . . ^ } and + 1 . . . N}, 
respectively, see Fig. [T] The coupling is given by a hop- 
ping term t' H between the left- and right-hand sides of 
the system (i.e., between site N/2 and site N/2+ 1) and 
an additional coupling V H ( — l/2)( fijv/2+i — 1/2). 
Here we discuss the homogeneous system only (t' H = in 
and ^h = Vh) but our approach can be easily extended 
to systems with a weak link t' H < tn representing the 
tunneling through a nanostructure2£ or with a site rep- 
resenting a quantum dot, such as the Interacting Reso- 
nant Level Model (IRLM)^, as well as to systems with 



a few additional sites intercalated between both leads 
and representing nanostructures with internal degrees of 
freedom^ 

The current operator between the pair of sites (fc, fc+1) 

is 



3k 



c k c k+l c k+l C k 



(5) 



For a given (time-dependent) quantum state we define 
the current flowing between both halves of the system 
(see Fig. [T]) as the expectation value of the current oper- 
ator for the site pair in the middle of the system 



J(t) = (3N/2) 



We note that 



(6) 



(7) 



where 



QUt) 



N/2 



and Q R (t) 



k=l 



N 

k=N/2+l 



(8) 

are the (time-dependent) charges in the left- and right- 
hand halves of the chain, respectively. As the number of 
particles n is conserved in our models, Ql{€) + Qn(t) = 
—e ■ n = const. The stationary current is a constantly 
flowing current in an infinitely large system after the set- 
tling time 



J 



lim lim J{t) 



(9) 



and according to the TLL theory in the linear regime 
(small Ae) it is given by£ 



h h 



for V > 0. 



(10) 



We set % = a L = 1 for all numerical simulations and, if 
units are not given explicitly, e = h = 1. 

In this work both setups used to generate a current 
(see the next section) lead to an overall half-filled sys- 
tem. Thus in the weak potential bias regime both sys- 
tem halves remain approximately half-filled at all times. 
Consequently, we expect that the current is given by the 
Luttinger liquid prediction (fTU| together with the Bethe 
Ansatz result ([3]) in the linear response regime. 



B. Setups for non-equilibrium simulations 

We employ two different setups to generate currents in 
the lattice models.— In the first one (I) we prepare the 
system at time t = in the ground state \(f>(Ae ^ 0)) of 
the Hamiltonian H = Hq + Hb (i.e. with potential bias), 
see Fig. [5] For later times t > we let the system evolve 
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FIG. 2. Setup (I): The charge reservoirs (two halves: source 
and drain) have different potentials but are coupled for t — 
while for t > the potential difference is set to Ae = 
instantaneously. 



according to the Hamiltonian Hq (i.e. without potential 
bias) 



\ip{t > 0)) = exp -i 



H t 



|</>(Ae^0)>. (11) 



This setup describes an inhomogeneous initial state with 
more particles in one half of the system than in the other 
one. Thus particles flow from one side to the other one 
for t > 0. It corresponds to a one-dimensional scattering 
experiment in which particles are emitted on one side of 
the system with energies between [—Ae/2, Ae/2], scat- 
tered at the junction between both system halves, and 
then (partially) transmitted to the opposite side. This 
picture of transport through junctions is often used in 
theoretical investigations. 

In the second setup (II) we prepare the system at time 
t = in the ground state |f6(Ae = 0)) of the Hamiltonian 
Hq (i.e., without potential bias). For later times t > 
the time evolution of the system is determined by the 
Hamiltonian H = Hq — Hb, i.e. with a potential bias 
that causes the current to flow in the same direction as 
in setup (I) 

m > 0)) = exp (-i {H ° ~ HB)t ) IflAe = 0)). (12) 



Setup (II) describes the evolution of an initially homoge- 
neous state under the influence of a potential gradient, 
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FIG. 3. Setup (II): The charge reservoirs (two halves: source 
and drain) are in equilibrium and coupled with tn and Vh for 
t = 0. For t > a potential difference Ae > is applied. 



see Fig. [3l Thus it corresponds more closely to the actual 
experimental situation with a voltage source generating 
a current in a conducting wire. 

We should point out that in setup (I) the two leads are 
decoupled with respect to the Coulomb- interaction for 
the calculation of the ground state. Our tests have re- 
vealed that otherwise a strong dependency of the station- 
ary current on the system size appears, that is, for smaller 
system sizes the stationary current becomes higher but 
for N — > oo it approaches the same constant value as 
the value one gets with = for t = 0. Choosing 
= for the computation of the ground state therefore 
decreases the finite-size error. 

Generally, currents (j6]) calculated with the states (fTTj) 
and (fT2|) are different. In the strong-bias limit |Ae| 3> 
in, Vjj it is easy to show that the steady-state current 
remains finite for the first setup while it vanishes for the 
second one. Recently, it has been reported that initial 
conditions (quenching an interaction term or a tunnel- 
ing term) can also alter the steady-state current flowing 
through a quantum point contact between two TLL leads 
which have been driven out of equilibrium by an exter- 
nal bias^i In the weak-bias limit |Ae| <C £hj however, a 
simple perturbation calculation shows that both setups 
yield the same linear response for the stationary current. 
Thus in this regime both setups can be used indifferently 
but in the non-linear regime we must distinguish them. 

In most theoretical studies the potential bias switch- 
ing is not instantaneous but adiabatic. This can also be 
used in numerical simulations, for instance, see Refs. 
and [1?]. Our tests do not reveal any significant differences 
for the steady-state current depending on the switching 
rate as long as the potential changes in a time scale which 
is much smaller than the time scale associated with the 
motion of particles from one reservoir to the other one. 
Since our extrapolation method is designed to work with 
an instantaneous change in the potential difference and 
as numerical simulations are simpler with it, we prefer 
this approach. 



C. One- particle equation of motion 

Without Coulomb- interaction Vh(= V^) = the 
spinless fcrmion model, described by the Hamiltonian 
((T|), reduces to the tight-binding model (non- interacting 
fermions without spin degree of freedom). A one- 
dimensional chain in the tight-binding model can be de- 
scribed by the single particle reduced density matrix 



g l3 (t) := (*(i)| C t C ,|*(t)> 



(13) 



where the time evolution is given by the one-particle 
equation of motion 



dt 



G(t) 



(14) 



-fft>o denotes the single particle Hamilton matrix of size 
N x N with which the system is evolved in time. More 
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precisely -ff t>0 i s t ne one-particle representation of Hq 
or H + Hb depending on the specific setup. The parti- 
cle number expectation values coincide with the diagonal 
terms of the reduced density matrix 

(n k (t)) = <*(t)|4 Cjfc |*(t)> = g kk (t) (15) 

and the expectation value for the current operator ([5]) 
from site k to k + 1 can be taken from off-diagonal entries 

&(*)} = i^[G k , k+ i(t) - Q k+ i, k {t)]. (16) 

Consequently, the dynamics of a tight-binding chain can 
be computed numerically without additional truncation 
error with a runtime of 0(N 3 ) for any point in time. 

D. TEBD method 

The TEBD metho d 18 ' 19 is based on a specific repre- 
sentation of quantum states by matrix product states 
(MPS). For an iV-site lattice it is generally written 

I*) =VrM j ' i A [1] r [ %' 2 A [2] ... 

_ . r [JV- A tAT-l] r [iV] J V| iii2 _ _ Jn) 

where \j1j2 ■ ■ ■ Jjv) designs the states of the occupation 
number basis, X^(k = 1, 2, . . . , N — 1) are positive def- 
inite diagonal matrices and F™ k (k = 1,2,..., TV) are 
matrices satisfying orthogonality conditions 

(rW jk V ^-^y T [k]jk = I, 

jk 

Y^T^(\^) 2 (Y^y =1. (18) 

jk 

Sums over an index j k run over a complete basis of 
the site k (for instance 2 states for the spinless fermion 
model). Every quantum state of the Fock space associ- 
ated with a finite lattice can be represented exactly in 
this form if the matrix dimensions can be as large as the 
square root of the Fock space dimension (for instance 
2 N / 2 for the spinless fermion model). In numerical com- 
putations, however, the matrix dimension must be kept 
smaller than a relatively small upper limit Xc- Fortu- 
nately, for many one-dimensional systems this trunca- 
tion is possible and can lead to a dramatic computational 
speedup while keeping the error in computed obscrvablcs 
conveniently low. This is done by using the Schmidt de- 
composition of the density matrix at each bond to calcu- 
late the matrices A^ from its eigenvalues and the T^ k 
from its eigenvectors. Given a bond k in the chain, the 
Schmidt decomposition at this bipartite split is defined 
by 

|*>= }T A^L 1 ^ 1 )!^ 1 -^ 1 }. (19) 

«fc=l 



The vectors arc the eigenvectors of p^" k \ the 

reduced density matrix of the left side of the split, the 
|^jfc+i..iv]^ correS p 0nc j m giy the eigenvectors of p^ k+1 -- N ^ 
for the right side, and the A^ fc the eigenvalues of both 
pi 1 -® and plfcH-N], with A Qfc > 0, ^* =1 A^ = 1. The 

Aq,,, are the matrix elements (A^ fc ]) Qfc and either of the 
|$ Qfc ) can be used to build the matrices T^ k from equa- 
tion (JT7J) . Thus, it is possible to keep only the largest Xc 
eigenvalues and throw away the rest while re-normalizing 
the state, such that the sum of the discarded values is 
smaller than an arbitrary error e 

E A L < e - (20) 

«t=Xc 

The best candidates are states for which the Schmidt 
dimension, and hence the dimension of the MPS matri- 
ces is low (like product states), or states for which the 
Schmidt coefficients display an exponential decay, where 
a large number of eigenvalues can be discarded without 
significant information loss. 

The great advantage of the TEBD algorithm is the 
possibility to compute the time-evolution of a state using 
a time-dependent Hamiltonian 

\*{t + St)) = e-* H ^ St \V(t)). (21) 

In a numerical implementation, St has to be discrete, 
such that the total simulation time r = n t ■ St, where n t 
is the number of time steps and St the numerical time 
step. The time-evolution can be used as well to calcu- 
late the ground state \ip gr ) of a Hamiltonian H . This is 
done by taking the time to be imaginary and projecting 
(effectively 'cooling') a starting state \ijjp) to the ground 
state of H 

e~ HT \ibp) 

ijj gr ) = fim ,, = \ VP .. . (22) 

Given a one-dimensional system of size N with nearest- 
neighbour interaction 

JV JV 

#" = E>f +E A t' +1] (23) 

1=1 1=1 

we can split the Hamiltonian into two sums = F + G 
over even and odd sites, with 

F = E (K[+K l 2 l+1 )= E F[l] > 

even I even I 

G = £(Al + /4< +1 )= £ GfW. (24) 

odd I odd I 

By using the Suzuki- Trotter decomposition for exponen- 
tial operators 

e ~i HSt = e -k(F+G)6t 

= e -iT«t e -iG5t e -i%8t + (St 3 ), (25) 
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FIG. 4. (Color online) Discarded weight D N / 2 for different 
parameters, solid black curve: Ae = 0.5tn,Vu = 0.8tn, 
dashed blue curve: Ae = 1£h,Vh = 1.6tH, dotted red curve: 
Ae = 24h,Vh = 0. The vertical lines indicate the runaway 
times described in the text. 



FIG. 5. Current for non- interacting fermions for setup (I). 
Shown are TEBD (dots on top of the dashed line) and exact 
results (dashed and solid lines) obtained with the equation of 
motion method for the single particle reduced density matrix 
(|13[) for different system sizes N. 



we can reduce the process of time-evolution to the suc- 
cessive application of operators which act only on two 
sites 

e ft 2 ot = e n 2 dt : 

even I 

e~i GSt = JJ e -* Gll]5t . (26) 

odd I 

The TEBD algorithm then proceeds to update the MPS 
representation by calculating the new Schmidt decompo- 
sition at the corresponding bond each time after a two- 
site operator is applied to it. This succession of two-site 
operators imposes a linear dependence of the computa- 
tional cost on the system size N. Each time a Schmidt 
decomposition is computed, we can truncate the dimen- 
sion of the MPS by keeping only a maximal number of 
eigenvalues, using a suitable limitation criterion. This 
is desirable, as the computational cost of both Schmidt 
decomposition and update of the MPS representation is 
0(xc 3 ) as explicitly shown in Refs. and [l9l The al- 
ternate updating of even and odd sites according to (|2"fJ)l 
makes the TEBD algorithm highly parallelizable using 
up to (N — l)/2 threads with a very low communication 
overhead. 

One obvious error source is the one stemming from 
the discretization of time in order to numerically com- 
pute equation (|21[) . so it is necessary to keep the actual 
time step St small enough to have a good approximation 
of H(t) in the interval [t,t + St]. But this error source 
is not inherent to the TEBD implementation. The main 
error sources in the algorithm are the Suzuki- Trotter ap- 
proximation and the Schmidt truncation. In order to 
improve the time step error, one can use higher dimen- 
sional Suzuki- Trotter formulas, at a computational cost 
which scales linearly with higher-order approximations^, 



or we can decrease the time step St, which also comes at 
a linear cost n t = ^. The dominating error in probably 
all setups of interest is thus the truncation error. It is 
also very difficult to compensate for this error, because of 
the 0(xc 3 ) scaling of the computational cost. As a trivial 
example, when starting a real time evolution with a state 
with Schmidt dimension smaller than Xc the simulation 
runs with a constant truncation error. As the simulation 
continues, there comes a point where more states than Xc 
are needed, and the truncation error quickly overcomes 
the Trotter error, which basically defines a runaway time 
for the simulation, as reported in Ref. H^. A reasonable 
and often used estimator of the truncation error is the 
discarded weight 

Q!fc=l 

where Dk denotes the discarded weight for a bipartite 
split at the £;-th bond and Xm < Xc is the number of kept 
eigenvalues. For our real time simulations, we use a max- 
imal Schmidt dimension in the range of 300 < Xc < 500, 
depending on the specific parameter combination. We 
use a site-dependent % c (site) and allow our simulations 
to adapt Xc(site) if test values (discarded weight, von- 
Neumann entropy) show the necessity to do so. We 
have to remark that for higher |Vh| and Ae the corre- 
lations (and thus the needed Schmidt dimension) within 
the spinless fcrmion model grow very quickly. Hence, the 
chosen Xe-range can lead to significantly larger errors in 
the regime |Vh| > 1.6£h and Ae > 3£h- 

Figure 2] shows the discarded weight for different pa- 
rameters for a split in the middle of the chain. One 
can see the runaway time indicated by the vertical lines. 
However, the maximal discarded weight over all simula- 
tions, sites and times up to 40h/tn we measured is yet 
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FIG. 6. Classical LC line with on-site applied potentials cj>i. 



smaller than 10~ 6 . Our simulations reveal that when 
taking a maximal Schmidt dimension of Xc = 600 in- 
stead of 300, the runaway time is shifted only about 20% 
to larger times and our quantities of interest are only 
slightly changed. Moreover, the stationary current com- 
puted with our method, which is described later, is only 
changed by less than 0.5%. 

Figure [5] shows the very good agreement between our 
TEBD simulations and the exact calculation using the 
single particle reduced density matrix (fF3")l for spinless 
fermions in the tight-binding model. We use a time step 
8t = 0.01 and a typical simulation runs approximately 
1/2 to 8 hours for a ground state calculation and about 4 
hours to (seldomly) 1 week for a real time evolution. We 
use a multi-threaded version of TEBD with an extremely 
low overhead (less than 1% for 10 processors and less 
than 5% for up to 96 processors). The overall cost of a 
TEBD simulation is slightly higher than that of a DMRG 
calculation on a single processor. If several processors 
are used in parallel the overhead of DMRG calculations 
become rapidly prohibitive, exceeding 100 % for as few 
as 4 processors^ Thus TEBD is already twice as fast 
as DMRG for 4 processors and the difference is likely to 
increase rapidly with the number of processors. We have 
used 12 threads for a TEBD simulation on average. 



III. CLASSICAL LC LINE 

In the following, our attention will be put on a classical 
model: the so-called LC line, shown in Fig. [5] The LC 
line has been under research as an electrical circuit im- 
plementation of a Toda chainiSr— , a nonlinear oscillator 
chain that, among others, describes the propagation of 
soliton waves. Further investigations of the LC line have 
so far been focused on the continuous and linear case4£ 
While this setup has therefore been analyzed for the non- 
linear and the linear and continuous case, our focus lies 
on a discrete and linear model for finite and infinite sizes. 

In the following, we show that the LC line serves well 
as the classical representation of a one-dimensional quan- 
tum wire described by the tight-binding Hamiltonian, 
while our main focus lies on the behaviour of the current, 
especially the influence of finite-size effects and station- 
ary values for the infinite-size limit. We present solutions 
for a setup in which an initial imbalance of charge car- 
riers on the condensators leads to an oscillating, rectan- 
gular current curve as the charge carriers move back and 



forth in the LC line: the formerly described setup (I). 
Scaling the system size N to infinity leads to an enlarge- 
ment of the square wave period, and for large times a 
stationary current / is flowing. Subsequently, a method 
is presented to compute the stationary current / for in- 
finitely large systems using finite-size simulation results. 
Based on the comparison of the classical and quantum 
mechanical model this approach will be applied to quan- 
tum systems. 

Figure IH] shows the classical, electronic setup of a one- 
dimensional wire that possesses many properties of a one- 
dimensional quantum chain. This LC line is a combina- 
tion of condensators and inductors whereas the present 
model is extended to enable on-site applied potentials fa . 
For this setup, the following relations can be derived us- 
ing elementary electrotechnical relations and Kirchhoff 's 
laws 

Ji = -CiUi ■ i = 1,2, ..,N 
Vi = Liii; i = l,2,..,JV-l (28) 
= Ii + J i+1 ; i = 0,l,..,JV-l 

and 

Vi = Ui-U i+1 +fa-fa +1 ; i = l,2,..,AT-l (29) 

where Iq = In = was defined. Vi in equation ([2"9")l 
results from the potential difference between the left 
(fa + Ui) and right (fa+i + Ui+i) connection of the i-ih in- 
ductor. Cj denotes the capacity of the i-th capacitor, Ui 
the voltage drop over the capacitor, Li the inductance of 
the z-th inductor and Ii and Ji are the currents as shown 
in Fig. |51 Qi is the charge of the i-th capacitor and the fa 
denote externally applied potentials. From the equations 
above follows 

I = -MI + $ (30) 

with 

and <f>i = (fa — fa + i)/Li. We consider an LC line with 
equal capacities and inductances 

L, = L and C t = C Vi (32) 



and temporally constant applied potentials $ = for 
which it follows 
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FIG. 7. Current l25o(t) for a classical wire of size N = 500 
through the inductor at position N/2. L = C = 1, Ql — 

Qo,i<N/2 = 1 and Qr = Qoj>n/2 = —1 



All values for the classical system ((J,-, t, L, C, etc.) are 
given in S.I. units unless otherwise specified. 

While the rectangular shape of the dominant oscilla- 
tion is understandable through the picture of charge car- 
riers moving back and forth in the line, the rapid oscilla- 
tion on top of the square wave is harder to explain and 
we will later take a closer look at it. The period T™ ax 
of the square wave can be calculated by using equation 
(f3"!)|) and PU]) and taking only the smallest frequency into 
account 



Timax 
cl 



sin i 



(41) 



For large system sizes (N ^> 1) this period is approxi- 
mately given by 



Timax 
cl 



(42) 



With Qo.i := Qi(t = 0) and 4>o.i := <f>i(t = 0) one has 

= jjj (Qo,i - Qo,*+i) + 2 (&M ~ 0<M+i) (35) 



and thus 



V2 



N-l r 



£ 



i 



(Qo,i — Qo,i+l) 

kin 



N 



(36) 



A. Stationary current 

In the following, we assume that initially (at t = 0) no 
current is flowing 



1,(0) =0 Vi 



Vfc 



(37) 



and that the chain is divided into two halves in which 
the charge carriers are uniformly distributed, i.e. 

Ql = <9o,i=i,..,f i Qr = Qo,i=f +i,..,jv 

and 4>i{t) = Q \/i,t, (38) 

which corresponds to the previously mentioned setup (I) . 

Figure [7] shows the current through the inductor Ljv/2 
according to the general formula for the current at site 
N/2 for any even N 



I N /2{t) = 



Ql-Qi 



N/2 

fe=i 



3in (vfc^) 



with 



/(2k - 1)tt 
V 2N 



(39) 



(40) 



and is thus linearly dependent on JV. In the second half 
period of the current in Fig.[7]one can see a beat upon the 
rapid oscillation. This is due to the unequally distributed 
charges which form an oscillating pattern and so create 
additional frequencies in the current. This effect also 
occurs in quantum wires, not only after the reflection 
at the right border but already from the beginning due 
to using an initial (ground) state with non-uniform local 
densities. 

In contrast to the oscillating current in a finite system 
where the charge carriers are scattered at the borders, 
we assume that the current in an infinitely large system 
is constantly flowing in one direction. To validate this 
assumption one has to execute the limit N — > oo first, 
and then take the limit t — > oo. In the limit N — > oo one 
gets from equation (|39| 



L 



>(*) = Jim lN/2(t) 




(43) 



with a = 2t/\[LC and ip = 2/\jLC . The solution is 
given by 



Ioo(t) 



(Ql - Qr) t 



[Jo {tpt) (2 - nHi (<pt)) 



2LC 

+ wJt {tpt) Ho ((ft)] 



(44) 



where J n (x) are the Bessel functions of the first kind and 
H^x) denote the Struve functions^ 

The current for the infinite system (|44[) is shown in 
Fig. [8] and it shows that for smaller times the curve for 
finite system sizes matches the one for the infinitely large 
LC line. Hence, a very accurate approximation of the 
value of a stationary current in an infinitely large system 
can be extracted from the value of a corresponding finite 
system result. This behaviour will be later used for the 
analysis of quantum systems. 
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FIG. 8. Solid curve: Current Iao{t) for a classical infinite 
wire through the inductor at the middle with L = C = 1, 
Ql = 1 and = — 1. The dashed curve is the current for 
a finite system with size N = 30. The constant line denotes 
the stationary value 1 against which the solid line converges. 



FIG. 9. Exact current Ioo(t) (solid line) and approximative 
current / a pp(i) (dashed line) with L = C = 1, Ql = 1 and 
Qr = -1. 



Using the approximations (|A1[) for the Bessel and 
Struve functions we find the value of the stationary cur- 
rent in the infinitely large system 

J cl = lim J eo (t)= Ql ~R r . (45) 
2y/LC 

The conductance G c i for the stationary case is given by 

with a local voltage drop 1/ = (Ql — Qr) jC in the mid- 
dle of the chain. 



of %i(x) in (|A2[) has an error of smaller order which 
presumably holds according to Ref. [45|. 

As shown in Fig. [5J expression (|4"T)l is a pursuasive 
approximation even for smaller times and thus a good 
explanatory basis for short-time effects in the LC line. 
It states a general 0(\/yfi) decay of the amplitudes of 
the rapid oscillation in the current curve which is the 
maximal order of the error for an approximation of the 
stationary current (|45l) using a finite-time (i.e. a finite- 
size) result. 

Instead of looking at the current through a single in- 
ductor, one can consider the current through two or more 
neighbouring ones. We will show that the short time be- 
haviour of this current has some advantages compared 
to the simple ijv/2(*) curve. We define the quantity of 
interest for even N as 



B. Finite-time and finite-size effects 

In order to describe finite-time effects in the infinite 
system we seek for a simpler current expression which 
does not only give the correct stationary value but also 
a good approximation for the short-time behaviour. Us- 
ing the asymptotic series expansions (|A2[) . ioo(i) can be 
transformed into the following form for i > 1 

- ^+^^^[(10* - 32)(sin 2 (v,t) 
+ sm(ift) cos(ipt)) 



+ 4(sin(</?t) - cos(pi))] + O 



(47) 



with ip = 2/\/LC. Idcv{t) denotes the deviation from 
the stationary current I c i from equation f|45|) . The order 
O (t^ 1 ) of the rest term is only valid if the approximation 



Ir*,N/2(t) = \l N /2(t) + \ (lN/2+l(t) + Jjy/a-l (*)) (48) 

which is equal to the expression \{lN/i{t) + -T/v/2+iW)) 
whereas for the latter expression the subsequently ex- 
plained effect also occurs for odd N for the substitution 
JV/2-> (JV-l)/2. In the limit N -> oo,i-> oo, I m)N/2 (t) 
converges against the same value as which implies 

that the stationary current is as well assessable via the 
current through two or more adjacent inductors. From 
(|33| one gets 



l m,N/2 



(*) 



h - Qr v^ 1 1 . , ,s . 2 f k7r 

— — > — sin (Witt) sin — 

NLC f^u k V k ' \2 



1 + cos 



kit 



(49) 



with expression (|36[) for bj- and the step distribution for 
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FIG. 10. Solid curve: Current 7 m ,oo(t) through two adja- 
cent inductors in the middle; equation (|51[) . Dashed curve: 
Current Iao(t) from equation (|44[) through a single inductor. 
L = C = 1, Ql = 1 and Q R = -1. 



the charges In the limit N — > oo it follows 



,(*) 



Sill 



sin (a;) 



cos 2 (a;)da;. 
(50) 

The integral can be solved in analogy to (|4"3")l with the 
substitution £ = sin (x) and gives 



.(*) 



^)(^ (^)--i,-i 



2, 
(51) 



This result, plotted in Fig. [TU] together with equation 
(HH for the current Ioo(t) through a single inductor, 
shows almost no rapid oscillations. An approximative 
expression for the rapid oscillations I r (t) m idev(i) can 
therefore be derived using the difference between the two 
current expressions ([51]) and (|44|) 



/r(*)=im,oo(*) -■?»(*) 



<3z, — Qr 



M<pt) 



t 



- }5 + <p)Ji((pt) 



(52) 



and using the approximations (|A2p for the remaining 
Bcssel functions one has 



-?r,app(i) 



(Ql - Qr) 



[(5(^) 2 + 8(^) 2 



32V^M) 5/2 t 

+ 15tpt - 6) COs(<^t) (5((/9t) 2 - 8(<y9t) 3 

- 15cpt - 6) sin(^)] . (53) 

Hence the rapid oscillations have a radial frequency ip = 
2/VLC = 27r/T c 7 n and thus 



T™ n = WW. 



(54) 



C. Determining the stationary current from finite 
system values 

There are various methods to calculate the stationary 
current from simulation results stemming from systems 
with finite size. The most obvious and easiest way to 
do so is to use the curve of a current through the mid- 
dle inductor (or the mean value of currents through two 
or more adjacent inductors) and to compute the mean 
value over a certain time interval in the quasi-stationary 
regime. In contrast to this approach, we will present 
a more analytical approach using the Fourier transfor- 
mation. Although one cannot prove that the following 
strategy can be applied to quantum mechanical systems 
with complicated interactions, several numerical and the- 
oretical indications are given in sections IIVI and [V] We 
start with a square wave described by its Fourier series 



AA ^ sin((2fc 

7T ^— ^ 



k=l 



2k- 1 



(55) 



where m determines the number of harmonics, A de- 
notes the value of the quasi-stationary plateau and x = 
27r/T^ lax where T]^ ax is the period of the square wave. 
This signal has certain similarities with the current curve 
In /2(f). It possesses a quasi-stationary regime which be- 
comes infinite for \ ~ ► (equals N — > 00 for the LC line) 
and shows a rapid and decaying oscillation on top of it. 
The main difference between the composition of the cur- 
rent curve (f3U)) and the square wave is the presence of 
sine-functions that enclose the (2k — l)-terms. 

The current in (|39|) is dependent on rjk and it is thus 
convenient to write I^/\(t) where % = sin( ^at^ )■ 
Choosing an r/k without the enclosing sine- functions gives 



r (2fc-i)» i 

rL 2JV J 

■*■ ?\r In 



L N/2 



(*) = 



N/2 

2(Ql-Qr) y; sin ((2fc -!)<//£) 
irx/LG ^ 2k -1 



(56) 



with ip' = ipTr/(2N) = ixj (Ny/LC). Since expression 
(|56|) now describes a square wave, one can compare it 
with (|55[) and get its quasi-stationary value 



A = 



Ql- 



= Ic 



(57) 



which is equal to the value that was derived for the sta- 
tionary current in equation (|45j). Remarkably the enclos- 
ing sine-functions do not change the stationary value, i.e. 
the expression (|56[) and the current in an LC line have 
the same stationary value for N — > 00. In case of ex- 
pression (|56j) this value is equal to the quasi-stationary 
value (f57|) which leads to the idea to calculate the quasi- 
stationary value from a given current signal in a finite 
system (using the procedure described below) and use 
it as approximation for the stationary current value for 
N -> 00. 

Applying the Fourier-transformation 



In/2{^) 



I N /2{t)e^dt 



(58) 
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term by term to expression (|39|) gives 



L N/2 



k=l 



, % . (Ql - Qr) S{lo + iprjk) - S(u - tprjk) 
(ui) = i ■ ) 



Vk 

(59) 

with equation (|40[) for r/k and ip = 2/V LC. If numeri- 
cal data with m values are given, one can only apply a 
discrete Fourier transformation which preserves the area 
underneath a delta peak, such that for any interval At 
containing a single delta function a ■ S(x — xq) the peak 
is of height a /At. Taking only the first (fc = 1) and 
only the positive delta-peak into account and using the 
discrete transformation 



m— 1 

f(j) = - E f(k)e~ lkjuJ0 with 



2tt 



UJQ 



k=0 



and j = 0, 1, m — 1, one gets for the current 



Ir 



Im[J w/ 2(wi)]JVAt / 7T \ 
; • sin I I 



\2NJ 



(60) 



(61) 



where Im[/jv/2(wi)] is the height of the first discrete peak 
(k = 1). In fact, one cannot expect to know the precise 
form of rjk for complicated quantum systems. Instead, 
one can take a first order approximative expression. For 
N ;§> 1 it holds sin (77^) ~ which implies that the 



•.2NJ 2N 

given curve is a square wave, and thus 



-Im[7jv/2(wi)]At 



with a relative error— 



\AL 



< 



8N 2 



(62) 



(63) 



It is important to remark that equation (|62p is not only 
an exact result for a square wave, but also an excellent 
approximation for the LC line current with maximal er- 
ror according to equation (|6"3")l . 

As our simulations reveal, the charge carriers in the LC 
line coming from the left half first reach the right border 
at t « T c 'f ax /4. After this time the finite system size has 
a significant influence on the simulation results. Apply- 
ing a discrete Fourier transformation means to assume 
that a periodic function is given. For that, and to min- 
imize the upcoming error for times larger than T™ ax /4 
the following procedure is proposed: From a given cur- 
rent curve lN/2(t) stemming from simulation results, first 
determine a position t = t s in the vicinity of T c I j lax /4 - 
but smaller than T™ ax / 4 - where the curve could be mir- 
rored along a parallel to the y-axis such that the loss of 
continuity is minimized, which is the case for instance at 
inflection points. Then consider a 4t s -periodic function, 
constructed as follows 



In(t,t s ) 



I N/2 (t) for < t < t s 

I N/2 (2t s -t) fort s <t<2t s 

-lN/2(t - 2t s ) for 2t s <t< 3t s 

-iN/z Ws-t) for 3t s < t < U s 



(64) 



< 




FIG. 11. Constructed signal Iiz(t,t s ) from simulation results 
together with the Fourier series representation of a square 
wave lZ m (t) from (|55p with m = 8 and an ideal square wave 
R(t). The Fourier transformation of Im(t,t a ) and expression 
(|62[) was used to calculate the amplitude of R(t). 



We then apply the Fourier transformation to the con- 
structed signal (|64|) and use equation (|62|) to calculate 
the corresponding stationary value of the current, sec 
Fig. QU 



D. A final note on setups 

Instead of initially dividing the system into two halves 
with different charges Ql and Qa, we can distribute the 
charge carriers homogeneously over the chain, apply a 
temporally constant voltage (</>o,l — </>o,k) & n d get the 
same result for the current as a function of time, which 
can be seen in equation (|36[) for b^. The reason for this 
lies in the linear character of our model which is as a 
consequence independent from the chosen setup. This is 
the main reason why we cannot fully compare the LC 
line with a quantum system, since preparing the classical 
system in setup (II) gives the same outcomes as for setup 
(I), in contrast to our quantum simulation results. 



IV. FINITE-SIZE EFFECTS AND STATIONARY 
CURRENT IN QUANTUM SYSTEMS 

A. Finite-size induced oscillation 

With both setups (I) and (II) one observes qualita- 
tively similar currents in the tight-binding model (the 
spinless fermion model without interaction Vqj = 0) as a 
function of time, as shown in Figs. 151 and IT2"1 First, there 
is a small transient regime for t < t a < 3h/tn with very 
rapid and dominant small oscillations. For long times 
t > tb ^> T max the current becomes very irregular be- 
cause of the progressive dephasing of moving particles. 
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B. Rapid oscillation on top of the square wave 



100 150 
Time t [h/t H ] 



FIG. 12. Current in the tight-binding model for a long time 
scale, calculated using the one-particle equation of motion. 
The initial state was prepared with a completely filled left 
half, a completely empty right half and tn — for all bonds. 



Between t a and tb we observe an approximately rectangu- 
lar wave with a period 7 imax which diverges with increas- 
ing system length (the corresponding leading frequency 
o; max = 27r/T max converges towards in the thermody- 
namic limit), see Fig. [12j The period of the rectangular 
oscillation is given by (Appendix [B]) 



T max _ h . _ for N > L 



(65) 



For spinless fermions, using a semi-classical picture, 
particles first flow from one side of the system to the 
opposite side because of the inhomogeneous density (first 
setup) or the potential difference (second setup). Then 
they are reflected by the hard wall represented by the 
chain edge. As there is no dissipation in our model, all 
reflected particles flow back with the same velocity in the 
opposite direction until they reach the other chain edge 
and are again reflected and so on. 

The progressive degradation of the rectangular signal 
can be understood using the same picture. First, all par- 
ticles flow in the same direction but, as they have differ- 
ent velocities, they progressively come out of phase. For 
long times t 3> tb, which can be checked up to the numer- 
ical double limit t « 10 300 for non-interacting fermions 
(Vh = 0), our simulations show that the current does 
not go to zero but continues to oscillate with a period 
T max . This can be understood in the picture of the clas- 
sical LC line as well as for spinless fermions through the 
dominant amplitude of the o; max oscillation. It is basi- 
cally the same effect as if all contributions to a series 
representation of a square wave (|55[) would 'randomly' 
come out of phase. In that case, the contribution with 
the lowest frequency o; max , which is the frequency of the 
rectangular oscillation, determines the frequency of the 
remaining irregularly shaped curve. 



Setting the two expressions (|4lj) and ([65]) for the peri- 
ods of the classical and quantum mechanical system equal 
one gets 



for N > 1. 



(66) 



Further comparison with equation (|54j) gives for the pe- 
riod of the rapid oscillation 

T min = ny/LC « — for TV > 1. (67) 
4i H 

Instead of preparing the system in a ground state of 
a Hamiltonian with applied potentials, one can initially 
decouple all sites and fill the left half of the system with 
(uncorrelated) particles, leaving the right half empty. 
This almost corresponds to choosing Ae = 4in, insofar 
as a ground state with this applied voltage and coupled 
reservoirs still has some correlations and unequally dis- 
tributed on-site particle densities in it. The resulting 
current curve in Fig. [T^] in the first half period closely 
resembles the curve shown in Fig. [7] for the LC line. 

The magnification in Fig. [12] shows a rapid oscillation 
with a period according to equation (|6"7) . in contrast to 
the oscillation shown in Fig.[T3]for setup (I) which is 'dis- 
turbed' by oscillations stemming from correlations and an 
unequal distribution of particles in the ground state due 
to a lower applied voltage Ae < 4tn- Although the oscil- 
lation in Fig.H^for setup (II) is very regular, it does not 
stem from the same 'classical' origin, which is confirmed 
by the fact that on the one hand two adjacent currents 
do not cancel out and on the other hand the oscillation 
does not fulfill equation ([57]). 




20 

Time t [H/t H ] 



40 



FIG. 13. Current in the tight-binding model for Ae = 2tu, 
calculated using the one-particle equation of motion. The 
dashed lines are currents through a single bond in the chain 
while the solid ones represent the mean value of the currents 
through two adjacent bonds in the middle. 
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FIG. 14. Solid line: Jf(£) from equation (|68[) . dashed line: 
corresponding approximation JF,a PP (t) from (1691) . 



Since in the following we are more interested in the 
stationary current flowing from an infinitely large source 
(the left half of the chain) to a an infinitely large drain 
(the other half), we reference to Ref. [46| for a discus- 
sion of time-dependent currents and to Ref. HE] for a de- 
tailed analysis of short-time and finite-size effects (in the 
IRLM). 



C. Applicability of the extrapolation approach to 
quantum systems 

Although we have seen the close connection between 
the tight-binding system for initially uncorrelated parti- 
cles and the classical LC line, there exists a major dif- 
ference in the behaviour of the respective currents: the 
damping of the rapid oscillations. While for the classi- 
cal LC line a 0(t~i) dependence of the amplitudes ac- 
cording to equation (|47p is predicted, this decrement is 
different for a quantum system. 

The expectation value of the current for setup (I) and 
for initial conditions like in Fig. [T^] (one completely filled 
and one completely empty half) is for N — > oo given 
byiL2£ (Appendix P 



eta 



Mt) = [(Mut)) 2 + (JiM)) 5 



(68) 



with u> = 2tn/h. Using the asymptotic series expansions 
12[) yields for the current 



^F,app(i) 



eta 



16h(ujt) 2 
+ 4sin(2wi)] 



[5 - 32wtcos(2wf) 
4ei H 



(69) 



which is plotted together with expression ([68]) in Fig. [MJ 
The curve highly coincides with the one shown in 
Fig. [T2] for t < T max /4 which confirms the assump- 
tion that the borders in the finite quantum system for 



setup (1) only significantly change the current after a 
time T max /4, determined by the velocity of the parti- 
cle density front wave^ which moves from the middle 
of the chain to the borders. The quasi-stationary value 
J = Aetn/h from ([69]) suits equation ([10]) for Ae = 4i H 
and the amplitudes of the rapid oscillation are damped 
with 0(l/t). Moreover, expression (jbTj) for the period of 
the rapid oscillation - which was derived only by compar- 
ison with the classical model - is confirmed by equation 

Since the latter analysis shows a 0(l/t) damping of 
the rapid oscillation, the open question at hand is, if 
it is reasonable to use our extrapolation method for a 
quantum system. In fact, a pure square wave described 
by its Fourier series (|55[) has also a damping of the rapid 
oscillation of 0(l/t) and highly resembles the current 
described by ([6"5]) or depicted in Fig. IT2"1 when choosing an 
appropriate amplitude and number of harmonics. Since 
the mathematical error (|63[) of our method stems from 
the difference of the analyzed curve to a square wave, the 
extrapolation approach provides an even smaller error 
than (|63p for a quantum system. 



V. SPINLESS FERMION MODEL 

In this section we discuss our results for the station- 
ary current and the finite-system period in the spinlcss 
fcrmion model. The non-linear dynamics of this model 
and the closely related spin i XXZ chain have been 
studied in several works using td-DMRG mcthod a 23 i 39 ' 49 . 
Nevertheless, the steady-state transport properties of the 
spinlcss fcrmion model have not been investigated yet. 



A. Current-voltage characteristics 

Firstly, we remark that the value of the potential bias is 
given by the parameter Ae of the Hamiltonian from equa- 
tion (U), and it is not the result of a quantum measure- 
ment of an observable. Thus it corresponds to an external 
field applied to the system, and in an interacting system 
it is not always equal to the actual potential difference 
which is measured in experiments Moreover, in ex- 
periments the sample (quantum wire) must be connected 
by leads to the voltage source and measurement appa- 
ratus. This modifies the low-energy, low-tcmpcraturc 
behaviour of the sample below some crossover energy 
scale^^ Concerning this matter, the linear conductance 
of a TLL attached to one-dimensional leads has been cal- 
culated in several works£i~— It has been shown that the 
conductance of this setup is given by e 2 /h (which cor- 
responds to equation (|10[) with K = 1), independently 
of the Coulomb interaction Vh within the wire, and that 
this result coincides with experimental outcomes. 

Our analysis of on-site particle densities shows that 
the densities for the first and last site of the chain only 
change after T max /4. This is a major reason why we only 
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FIG. 15. (Color online) Current-voltage curve of the spinless 
fermion model with setup (I) for several positive Vh- For 
Vh = the current has been computed using TEBD and the 
one-particle equation of motion (eom). The lines are guides 
for the eyes. 



take values of the current for times smaller than T max /4 
into account, as described in our extrapolation approach. 
To test the overall procedure for computing the station- 
ary current we have first calculated the current-voltage 
characteristic of a non- interacting chain (Vh = 0). This 
characteristic is known exactly for the first setup 47 ' 48 : 
J = {e/h)Ae for \eV\ < 4i H and J = (e/h)4i H for 
\eV\ > 4iu, with Ae = \eV\ where V is the voltage bias. 
For the second setup, exact results were calculated nu- 
merically using the one-particle equation of motion with 
a system size of N = 1000. 

We have found that our procedure yields stationary 
current values which agree with an overall error of less 
than 5% with the exact results. As a second test we 
have calculated the current-voltage characteristic of the 
IRLM. The generic shape of this curve is known from a 
field-theoretical analysis, and highly accurate td-DMRG 
results are available for a quantitative comparison^ In 
that regard, we have found that the results from our 
method agree very well with the numerical outcomes 
from Ref. [l2|- We nevertheless have to remark that we 
found the IRLM far easier to simulate than the spinless 
fermion model, since for higher |Vh| the correlations (and 
thus the needed Schmidt dimension) within the spinless 
fermion model grow more quickly than in the IRLM. As 
already mentioned in section III Dl the parameters used 
for our simulations can result in larger errors for high 
voltages and Coulomb interactions in the TEBD data for 
the current as a function of time. Additionally, finite- 
size effects also grow with |Vh| and Ae. This can lead to 
less accurate extrapolations for the steady-state current 
in that regime than for non-interacting fermions and the 
IRLM. 

We now discuss the current-voltage characteristics 
which have been obtained with the methods described 
above. In setup (I) the current of the non-interacting 
system (Vh = 0) is strictly proportional to Ae = \eV\ up 
to the band width 4£h and then remains constant] 47 ' 48 



FIG. 16. (Color online) Magnification of figure [151 The lines 
indicate the linear response according to the TLL theory (|10|) 
for small Ae. 



For Vh > we see in Fig. [TS] that the current increases 
sub linearly with Ae for a fixed interaction strength Vh 
and that it decreases monotonically with increasing Vh 
for a fixed potential bias Ae. The magnification in Fig. 1161 
shows a comparison with the linear response in the TLL 
theory (fTQj) for small Ae, using the Bethe Ansatz solution 
parameters (J3J). The good agreement confirms the valid- 
ity of our approach. Obviously, an increasing Vh does not 
only reduce the current but also the range of the poten- 
tial bias Ae for which the linear response approximation 
is accurate. 

In Fig. [17] we observe a behaviour for attractive inter- 
actions Vh < which is similar to the non-interacting 
case. The current increases almost linearly with Ae up 
to some Vn-dependent threshold potential where it sat- 
urates at its maximal value. Increasing |Vh| scales the 
maximal current down but causes a higher linear con- 
ductance for small Ae. Figure [TBI shows that for negative 
interactions there is also an excellent agreement between 
our results and the linear response in the TLL theory 
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FIG. 17. (Color online) Current-voltage curve of the spinless 
fermion model with setup (I) for several negative Vh. The 
dashed lines indicate the theoretical beginning of the current 
plateaus according to (|70|l . The solid lines are guides for the 
eyes. 
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FIG. 18. (Color online) Magnification of Fig. El The lines 
indicate the linear response according to the TLL theory (|10|) 
for small Ae. 



(ITO)) for small Ae. Again increasing |Vh| seems to reduce 
the range of the potential bias Ae for which the linear 
response approximation is accurate. 

The occurrence of current plateaus for large potential 
differences can be easily understood. In setup (I) one 
half of the chain initially contains more particles than the 
other one due to the applied potential bias Ae. When this 
bias is large enough the initial state consists of one com- 
pletely filled and one completely empty half and it does 
not change for higher Ae. In the non- interacting case this 
saturation occurs exactly at Ae = 4£h' 47 ' 48 When we add 
an attractive interaction Vh < the particles are more 
likely to stick together and therefore even more particles 
gather in one half for the same applied potential. Thus, 
saturation (i.e., an initial state with one completely filled 
and one completely empty half) is reached for smaller Ae 
as |Vh| increases as seen in Fig.[T7J We do not understand 
yet why the plateau heights (i.e. the maximal current) 
are lowered by increasing |Vh|- For a repulsive interac- 
tion Vh the effect is reversed and less particles gather in 
one half of the chain for a given potential difference when 
Vh increases. Thus we expect that saturation occurs at 
higher values Ae > 4£h beyond the potential range shown 
in Fig. 1151 We can easily determine for which potential 
difference Ae saturation occurs. Removing a single parti- 
cle from the completely filled reservoir or adding one par- 
ticle to the empty half costs an energy Ae/2 — Vh — . 
Thus saturation occurs if 



Ae > 2V" H +4t H . 



(70) 



Figure [17] shows that this approximation fits well to the 
numercial data for Vh < whereas the saturated regime 
according to ([70]) for Vh > lies outside the potential 
range shown in Fig. 1151 as mentioned above. 

While for setup (I) the differential conductance is al- 
ways positive (or zero after saturation), in setup (II) we 
observe a negative differential conductance. Figure [19] 
shows our results for setup (II) with Vh > 0. The cur- 
rent seems to vanish for very large potential biases as 
predicted by strong-coupling perturbation theory for this 



FIG. 19. (Color online) Current-voltage curve of the spinless 
fermion model with setup (II) for several positive Vh- For 
Vh = the current has been computed using TEBD and the 
one-particle equation of motion (eom). The lines are guides 
for the eyes. 



setup. We note that the current becomes systematically 
weaker with increasing Vh for a fixed small potential bias. 
For larger Ae the behaviour of I as a function of Vh is 
no longer monotonic. In addition, we see that the onset 
of the negative differential conductance (i.e., the position 
of the maximum of I as a function of Ae) shifts to higher 
values with increasing Vh- The magnification for small 
Ae, Fig. [201 confirms again that our results agree with 
the TLL theory in the linear regime and that the range 
of the linear response regime shrinks with increasing in- 
teraction strength. 

Results for attractive interactions Vh < are shown 
in Fig. [2U There we clearly observe that the current 
vanishes for large potential differences as predicted by 
strong-coupling perturbation theory for setup (II). We 
also see that the current becomes systematically larger 
with increasing interaction strength |Vh| for a fixed small 
potential bias. For intermediate Ae the behaviour of I as 
a function of Vh is not monotonic and for large enough 
Ae > 2£h the current decreases systematically with in- 
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FIG. 20. (Color online) Magnification of Fig. \M The lines 
indicate the linear response according to the TLL theory (|10|) 
for small Ae. 
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FIG. 21. (Color online) Current-voltage curve of the spinless 
fermion model with setup (II) for several negative Vh- The 
lines are guides for the eyes. 



creasing interaction strength. Moreover, we see that the 
onset of the negative differential conductance occurs at 
lower values of Ae for increasing |Vh|- Finally, Fig. |2"21 
shows that in this case our results also agree with the 
TLL theory in the linear regime, but it seems that the 
interaction reduces rapidly the range of potential biases 
for which the linear response is accurate. 

The negative differential conductance which is ob- 
served in setup (II) for all Vh can be understood in terms 
of the overlap of the energy bands for elementary excita- 
tions. Our simulations show that there is no significant 
difference in the current curve whether one initially does 
or does not couple source and drain (i.e., the two halves 
of the system). For the following considerations we can 
therefore regard source and drain as initially decoupled, 
having two separate but equal energy bands. Then we 
apply a step-like voltage in setup (II) and thus shift the 
two bands against each other. For Vh = only over- 
lapping occupied one-particle states in the source and 
unoccupied ones in the drain contribute to the current. 
While this overlap rises with increasing |Ae| from to 
2£h, it diminishes from |Ae| = 2in and reaches zero at 
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FIG. 23. (Color online) Period T max of the rectangular cur- 
rent oscillation in the finite spinless fermion model as a func- 
tion of the inverse renormalized Fermi velocity v~ from ex- 
pression ([2]) . The dashed lines are described by equation (|71|) • 



4<h- As a result, the current is maximal for Ae « 2iu and 
vanishes for large Ae. Similarly, we can understand the 
non-monotonic behaviour in interacting cases (Vh 0) 
in terms of the overlap of the elementary excitation bands 
in the spinless Luttinger liquids in the two halves of the 
system. However, the effective bandwidth is renormal- 
ized exactly as the Fermi velocity in equation There- 
fore, the maximum of the current is reached for a smaller 
potential bias Ae when Vh becomes negative as shown 
in Fig. EH and shifted to a higher potential bias when 
Vh increases as seen in Fig . [T51 Our conclusion agrees 
with the findings in Ref. [36J where it is shown that within 
a similar one-dimensional model a negative differential 
conductance is mainly caused by finite electrode band- 
widths. 



B. Influence of Vh on the period T r 
quantum system 



in the finite 
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FIG. 22. (Color online) Magnification of Fig. ED The lines 
indicate the linear response according to the TLL theory (|10[) 
for small Ae. 



Our results show a further effect of Vh on the finite- 
system current. While for Vh > the period of the rect- 
angular oscillation becomes smaller, it grows for negative 
interactions. 

Figure [23] shows the period T max of the rectangular 
current oscillation in the finite system as a function of the 
inverse renormalized Fermi velocity w _1 from expression 
@ • The T max values are mean values from various Ae for 
each setup. These values agree perfectly with the dashed 
lines in Fig. [231 which are given by the equation 



T max (V H ) 



^(Vh = 0) 
v(V H ) 



T max (V H = 0). 



(71) 



Thus, the period is fully determined by the time-scale of 
the non-interacting system T max (VH = 0) from (|B4[) and 
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the renormalizcd Fermi velocity v from expression p]). 



VI. SUMMARY AND CONCLUSION 

We have presented a method for investigating the 
non-linear steady-state transport properties of one- 
dimensional correlated quantum systems. First, we have 
analyzed finite-size effects and transient currents in a re- 
lated classical model, the so-called LC line. This ex- 
act analysis shows how steady-state currents in the ther- 
modynamic limit can be calculated from transient cur- 
rents in finite systems. We have found that for currents 
in quantum systems, in particular the spinless fermion 
model, finite-size effects can be understood from the be- 
haviour of the LC line. 

On the basis of the strong connection between the 
finite-size behaviour in classical and quantum models we 
have developed a method to extract a stationary cur- 
rent from simulations of finite-size quantum systems. We 
only need numerical data I(t) up to a time of the or- 
der of T max /4 where T max cx N is the period of the ap- 
proximately rectangular signal I(t) in a quantum system 
of size N. The non-equilibrium dynamics of correlated 
quantum systems has been calculated using the time- 
evolving block decimation method (TEBD). We have im- 
plemented a multi-threaded version of TEBD with an ex- 
tremely low overhead (less than 1% for 10 threads) which 
allows us to simulate the non-equilibrium dynamics up to 
the time-scale 2T max . 

We have determined the full I-V characteristic of the 
spinless fermion model with nearest-neighbour hopping 
tu and interaction Vh using two different setups to gen- 
erate currents. In setup (I) the initial state has different 
particle numbers in its two halves due to an applied po- 
tential difference while the system evolves in time with 
an overall equal on-site potential. In setup (II) we cal- 
culate the initial state without a potential difference but 
turn it on for the real time evolution. For non-interacting 
fcrmions (Vh =0) our outcomes agree with the analytical 
solutions and with the results of the one-particle equa- 
tion of motion method. Additionally, we have checked 
that our results coincide with the field-theoretical analy- 
sis and td-DMRG simulations^ for the IRLM. 

For interacting fermions we have found that the steady- 
state current is independent from the setup in the linear 
regime \eV\ -C 4£h- Moreover, our numerical data agree 
with the predictions of the Luttingcr liquid theory com- 
bined with the Bethe Ansatz solution. For larger poten- 
tials V the steady-state current depends on the current- 
generating setup. This difference is well understood for 
non-interacting system o 35 ! 36 but can only be explained 
qualitatively for our interacting system. With setup (I) 
we have found that the current increases with V and sat- 
urates at a finite value when the potential difference is 
high enough to separate the initial state in one filled and 
one empty half. With setup (II) we observe a negative 
differential conductance that can be understood in terms 



of the overlap of the elementary excitation bands in the 
spinless Luttinger liquids in the two halves of the system. 
Both effects - the plateaus and the negative differential 
conductance at large potential bias - are due to the finite 
bandwidth of the system which agrees with the findings 
in Ref . [36l. However, in our case the interaction renormal- 
izes the effective bandwidth exactly as the Fermi velocity 
in equation ^ so that the current maxima and the onset 
of the plateaus depend on the strength of the interaction 
Vh- The change in the typical time-scale T max also solely 
depends on the renormalized Fermi velocity. 

The methods presented in this work can be applied to 
other systems, such as models of quantum dots or wires 
coupled to leads including electronic and also bosonic de- 
grees of freedom. We have already checked the applicabil- 
ity of these methods to the IRLM, the one-dimensional 
Hubbard model away from half-filling and also two-leg 
ladder systems. Hence, we believe that they will be very 
useful to study non-linear transport properties of cor- 
related low-dimensional conductors. Although we have 
mostly tested our appraoch in the limit of transparent 
coupling between source and drain, we think that it can 
also be applied when the hybridization between leads is 
very small. This issue will be examined in future works. 
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Appendix A: Approximations of the Bessel und 
Struve functions 



The Bessel functions of the first kind and Struve func- 
tions have the following approximations which become 
exact for x — > oo 



v^ OT (*-i) + °G)' 

— COS (x - -7T ] + 0[ - 

7TX \ 4 ) \X 



2 / 7T 
SIB \ X 

7TX V 4 



O 



Jo(x) 
Ji{x) 
Ho{x) 



Using these expressions together with equation (|44|) 
yields the stationary current value (|45[) but does not pro- 
vide an approximate expression for the short-time be- 
haviour of the current, since all time-dependent terms 
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cancel out. Instead, the asymptotic series expansions 
from Refs. and [53 and the approximation of 'Hi(x)^- 



Jo(x) 
Ji(x) 
Ho(x) 

Hi{x) 



(8x - 1) cos(x) + (8x + 1) sin(a;) 

8V TTX 3 

(3 — 8x) cos(x) + (3 + 8x) s'm(x) 



I A, 

X2 
1 
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result in the expression f|47|) . 
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Appendix C: Current in the infinite one-dimensional 
tight-binding model 



The expectation value of the current for setup (I) and 
for initial conditions like in Fig. [12] (one completely filled 
and one completely empty half) is for N — > oo given 
byiLiS 



2et B 



J2 Jt(ut)Ji + i{ost) (CI) 



where k denotes the site, u = 2tn/h and Ji(z) are the 
Bessel functions of the first kind. Reformulating the sum 
and utilizing the Bessel recursion relation 



2/ 



Jl(z) = Ji+i(z) + Ji-i{z) 



(C2) 



Appendix B: Period of the rectangular oscillation in 
the tight-binding model 



gives for z ^ and k = N/2 



The time evolution of the single particle reduced den- 
sity matrix (| 14[) is given by 



Gijit) = cxp (-^e,ijftfc 9 (0)exp (Uktj (Bl) 

in the eigenbasis of a time-constant single particle Hamil- 
tonian H^K ^k.i denotes the i-th component of the fc-th 
eigenvector of H^ 1 ' and et the corresponding eigenvalue. 
A Fourier transformation gives 



Q kq {io) = S [e k -Sq- hu] Q kq {0). 



(B2) 



For the tight-binding model with zero on-site potentials 
one has 
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The period of the largest (rectangular) oscillation is then 
given by 



nil 



in sin 



N+l 
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%■ — for N > 1. (B4) 
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for Rc(^) > — 1 and z ^ w the current is given by 
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Applying l'Hospital's rule one gets 
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